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Abstract. An abstract should be given 



Monte Carlo techniques based on indivisible energy 
packets are described for computing light curves and spec- 
tra for 3-D supernovae. The radiative transfer is time- 
dependent and includes all effects of 0(v/c). Monte Carlo 
quantization is achieved by discretizing the initial distri- 
bution of 56 Ni into Af radioactive pellets. Each pellet de- 
cays with the emission of a single energy packet comprising 
7-ray photons representing one line from either the 56 Ni 
or the 56 Co decay spectrum. Subsequently, these energy 
packets propagate through the homologously-expanding 
ejecta with appropriate changes in the nature of their con- 
tained energy as they undergo Compton scatterings and 
pure absorptions. 

The 3-D code is tested by applying it to a spherically- 
symmetric SN in which the transfer of optical radiation is 
treated with a grey absorption coefficient. This 1-D prob- 
lem is separately solved using Castor's co-moving frame 
moment equations. Satisfactory agreement is obtained. 

The Monte Carlo code is a platform onto which more 
advanced treatments of the interactions of matter and ra- 
diation can be added. Some of these have already been 
developed and tested in previous papers and are summa- 
rized here. 
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1. Introduction 



radiative transfer 



In earlier papers, Monte Carlo (MC) methods were used to 
construct spectral synthesis codes for SNe in their photo- 
spheric phases. Because the emphasis was on coarse anal- 
yses of the spectra of newly-discovered SNe, major sim- 
plifying approximations were made to minimize comput- 
ing time and to ensure robustness. Thus, in the first code 
(Lucy 1987; Mazzali & Lucy 1993), the SN's atmosphere 
is spherically-symmetric and homologously expanding, the 
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radiation field is stationary (c = oo), continuum formation 
is confined to a sharply-defined lower boundary (Schuster- 
Schwarzschild approximation), and line formation results 
from the coherent scattering of this continuum as it prop- 
agates through the outer layers. Moreover, the stratifica- 
tions of temperature, ionization and excitation are derived 
not from first principles but from formulae that approx- 
imate the effects of dilution in the SN's extended atmo- 
sphere. 

In a subsequent paper (Lucy 1999b), two innovations 
significantly improved this code. First, coherent scattering 
was replaced by downward branching as the mechanism of 
line formation. Secondly, a computed spectrum's sampling 
errors were greatly reduced by using the formal integral 
for the emergent intensity instead of simply binning the 
escaping photon packets. 

Although these codes have proved their worth diagnos- 
tically, they cannot compute spectra for explosion models, 
a challenge that must be faced if explosion mechanisms 
and progenitor scenarios are to be confronted with ob- 
served spectra and light curves. Evidently, simplifying as- 
sumptions must be reconsidered in planning a more pow- 
erful and versatile code. 

First, the assumption of stationarity must be aban- 
doned. As Arnett (1980, 1982) long ago demonstrated for 
SNe both of types I and II, the luminosity L(t) at elapsed 
time t does not in general closely approximate the instan- 
taneous energy deposition rate by 7-rays emitted in ra- 
dioactive decays. Accordingly, the time-dependent diffu- 
sion of radiant energy through the expanding ejecta must 
be treated explicitly if we wish to compute the emergent 
luminosity density L v as a function of time. 

Similarly, artificially creating radiant energy by means 
of a continuum-emitting lower boundary is no longer ac- 
ceptable. The computational domain must be the com- 
plete configuration, not just a 'reversing layer'; and con- 
tinuum formation must be treated explicitly. Escaping ra- 
diant energy then derives ultimately from radioactive de- 
cays or from the energy content of the ejecta at ti, the 
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time at which the output from an explosion calculation is 
used to initiate the spectral synthesis code. 

Removal of the constraint of spherical symmetry is also 
highly desirable in view of accumulating observational ev- 
idence and theoretical arguments implying that most and 
perhaps all SN explosions are significantly asphcrical - see 
Wheeler (2004) for a recent review. 

Of the assumptions in the earlier codes, the only one 
retained is that of homologous expansion. This requires 
that the explosion calculation, which of necessity includes 
gas dynamics, be continued until all mass elements are to 
a good approximation coasting ballistically. 

With the assumptions to be dropped identified, the 
problem is now defined: to solve the time-dependent, 3- 
D NLTE transfer problem for UVOIR radiation in the 
homologously expanding ejecta of a SN, given the distri- 
bution of mass and composition at an initial time t\. Of 
course, this problem is coupled to a corresponding 3-D 
transport problem for the 7-rays emitted by radioactive 
isotopes. This coupling occurs via the equations describ- 
ing the deposition and degradation of 7-ray energy. 

In view of the magnitude of this problem, this paper is 
restricted to describing an exploratory MC code in which 
the UVOIR radiation's interaction with matter is governed 
by a grey absorption coefficient. This problem was previ- 
ously treated by Pinto & Eastman (2000) in order to ex- 
plore the sensitivity of type la light curves to parameters. 
Here the aim is to create a software platform onto which 
more realistic physics can be subsequently added and with 
which numerical techniques can be tested. 

Monte Carlo methods are a natural choice for this 
problem given the experience with the earlier codes. 
Moreover, MC methods are commonly favoured for trans- 
port phenomena in geometrically-complex configurations 
with no symmetries. Nevertheless, the conventional ap- 
proach in which the derivatives in the time-dependent 
transfer equation (RTE) are approximated by differences 
has already been partially implemented for 3-D SNe by 
Hoflich (2003). With regard to the application of MC 
methods, the most advanced work seems to be that of 
Kasen et al. (2004) who carry out 2-D timc-indcpcndcnt 
calculations to compute the observable characteristics of 
a type la model with a conical hole. 

2. Monte Carlo techniques 

In this section, after some preliminaries, recent develop- 
ments in MC technique that are relevant for a general 
spectral synthesis code are summarized. These remarks 
are not restricted to the grey case considered in this pa- 
per. 

2.1. Random numbers 

Random numbers are obtained with a double precision 
version of the routine ran2 of Press et al. (1992). Such 
numbers are always denoted by z, with each z denoting 
an independent call to ran2. 



2.2. Energy packets 

As in earlier codes, the MC quanta are energy packets. 
In general, these are referred to as ^-packets except when 
specifying the nature of the contained energy (Paper I, 
Sect. 2). Thus r-packcts and 7-packets are monochromatic 
photon packets of UVOIR radiation and of 7-rays, respec- 
tively. On the other hand, fc-packets and i-packets con- 
tain not radiation but thermal kinetic energy and ioniza- 
tion energy, respectively. In addition, a full treatment of 
radioactive decay and 7-ray deposition will require con- 
sideration of e~- and e + -packets containing non-thermal 
electrons and positrons, respectively. 

2.3. Discretization 

The SN is enclosed in a Cartesian grid containing I 3 iden- 
tical cubic cells. This grid expands with the ejecta so that 
each cell contains a fixed parcel of matter. Physical vari- 
ables do not vary spatially within cells. 

The expansion itself is discretized into N time steps. 
Expansion occurs in instantaneous jumps at times t n sep- 
arated by a constant value of A log t. During every time 
interval (t n ,t n+ i), the density p of each cell is held fixed 
at the value predicted by homologous expansion at time 

tn+l/2 = y/tn tn+1- 

Here, and throughout the paper, t is the elapsed time 
since explosion and is measured in the SN's centre-of-mass 
rest frame (rf). 

With this time-stepping, the calculation is divided into 
N separate MC experiments, each of which updates the 
radiation field throughout the grid from time t n to t n+1 
for n = 1,2,..., N. 

2.4. Relativistic terms 

Given the modest accuracy with which we can currently 
treat many of the physical processes in SN ejecta, all terms 
of 0(v/c) in radiative transfer could be neglected except 
for Doppler frequency shifts, which are essential in line 
transfer. Nevertheless, looking forward to a time when the 
relevant cross sections will be known to high precision, 
terms of 0(v/c) are now included. 

The basic transfer calculations are carried out in the 
rf, but transformations to and from the local co-moving 
frame (cmf) are convenient in treating interactions with 
matter. The motion of an £ -packet is thus given by r(t), its 
position as a function of time t in the rf Cartesian coordi- 
nate system with origin at the SN's centre of mass. When 
the grid attached to the ejecta expands instantaneously at 
times t n , the coordinates r and direction vectors fi of the 
^-packets remain unchanged. 

2.5. Energy conservation 

In the context of this discretization of space and time, 
energy conservation implies I 3 N constraints that the so- 
lution must satisfy. Thus, for each cell and every time step, 
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the net emissivity of radiant energy per unit mass must be 
balanced by energy created within the mass element after 
subtracting the increase dU in its internal energy and the 
pdV work done by gas pressure. 

With conventional transfer techniques, a solution sat- 
isfying this huge number of constraints - perhaps ~ 10 8 - 
has to be achieved iteratively - and convergence is never 
rapid. In contrast, if the quanta in a MC calculation are 
indestructible and indivisible ^-packets, local energy con- 
servation is automatically and rigorously obeyed, even if 
physical variables have not converged to their final values 
(Lucy 1999a). 

To a good approximation, energy conservation can be 
simplified to the condition of thermal equilibrium by ne- 
glecting the gas terms dU + pdV (Arnett 1980; Pinto & 
Eastman 2000a). With this simplification, the net emis- 
sivity for UVOIR radiation in the cmf is equal to the rate 
of energy deposition by 7- rays. 

Thermal equilibrium is rigorously obeyed by the MC 
calculation because the effective cmf emissivity j e ^{u) 
implied by using indivisible ^-packets is subject to the 
integral constraint 

4tt J j' eff (v)dv = 4tt J k u j' u du + 'H 1 (1) 

where k v is the cmf absorption coefficient per unit volume, 
J' v is the cmf mean intensity, and 7i_, is the cmf 7-ray 
energy deposition rate per unit volume. In contrast, when 
the RTE is solved conventionally, the cmf emissivity j v 
is given by the fundamental NLTE formulae relating the 
emission of photons by bound-bound (b-b) and free-bound 
(f-b) transitions to the level populations and the electron 
temperature. Eq. (1) is then only satisfied asymptotically 
as the iterative procedure converges. 

Note that, in contrast to j v derived from first priniples, 
the MC emissivity j e ^(u) is not defined by mathematical 
formulae. Rather it is defined operationally by the rules 
governing the re-emission of absorbed ^-packets (Papers 
I, II). 

2.6. Statistical equilibrium 

In addition to the constraints demanded by thermal equi- 
librium, an even larger number of constraints are required 
by statistical equilibrium. During each time step and for 
each cell, the number of excitations of an atomic level 
must balance the number of de-excitations. With accu- 
rate models for the atoms in the multi-species plasma, the 
total number of levels might be ~ 10 3 . This is therefore 
the number of constraints per cell per time step implied 
by statistical equilibrium. Combining this with the ear- 
lier estimate of I 3 N, we see that ~ 10 11 constraints must 
be satisfied in the course of computing NLTE light curves 
and spectra of a 3-D SN. 

Indivisible ^-packets allow the constraint of statistical 
equilibrium to be incorporated automatically and rigor- 
ously, even when excitation and de-excitation rates com- 
puted from basic formulae do not balance. This is achieved 



with the macro-atom formalism (Paper I), according to 
which the rate of de-excitation of an excited state is not 
determined by its level population but by the rates at 
which radiative and collisional processes from all other 
levels populate the level in question, cither directly or 
through the internal transitions of the macro-atom. 

If the radiation field is computed from the RTE, errors 
in level populations, such as those present prior to conver- 
gence, result in the non-physical creation or destruction 
of radiant energy This problem is a consequence of the 
generality of the RTE equation with the NLTE emissivity 
j v . Because of this generality, the NLTE RTE could be 
used to follow the time-dependent relaxation to statistical 
equilibrium. In contrast, with the macro-atom formalism, 
spurious sources or sinks of radiant energy are eliminated. 
This is achieved by effectively solving a less general trans- 
fer equation, namely one that incorporates the constraint 
of statistical equilibrium and one, therefore, that cannot 
be used for a relaxation calculation. 

A consequence of computing de-excitation rates with 
the macro-atom machinery is that these rates and the 
corresponding emissivity j e ^{v) are insensitive to errors 
in the populations of the emitting levels (Paper I, Sect. 
6). In contrast, with the RTE, errors in level populations 
translate directly into errors in j v and therefore in the de- 
rived radiation field. Evidently, the macro-atom approach 
is more tolerant of departures from the NLTE solution 
than is the conventional approach. 

As thus far developed, the macro-atom approach im- 
poses the constraint of thermal equilibrium simultaneously 
with that of statistical equilibrium. Thus when a fc-packct 
is created (Paper II, Sect. 4.3.1), it is instantaneously elim- 
inated (Paper II, Sect. 4.3.2). Thus, there is no net energy 
transfer to or from the thermal pool. Moreover, this holds 
even if the electron temperature does not imply thermal 
equilibrium. Note also that if a fc-packet is created by the 
de-activation of a macro-atom of one atomic species, it 
may be eliminated by the activation of a macro-atom of 
another species. 

2. 7. Non-equilibrium 

For SNe ejecta at late times , the assumptions of ionization 
and thermal equilibrium break down (Fransson & Kozma 
1993). But this does not preclude the use of indivisible £- 
packets, since these are fundamentally a means of tracking 
energy. Departures from equilibrium can be modelled by 
allowing for the finite life time of a non-radiative £-packet 
before it converts back into an r-packet. Thus ab — f pro- 
cess creates a fc-packet or an i-packet, but neither should 
be immediately eliminated if the recombination time is 
not short compared to the elapsed time t. Moreover, while 
awaiting elimination, the fc-packet 's energy declines due to 
pdV work. 

If, during such non-equilibrium phases, statistical equi- 
librium remains a good approximation for the excited lev- 
els of each ion separately, then this constraint can be im- 
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posed by introducing macro-ions adapting the procedures 
of Paper I. 

The above remarks strongly suggest that non- 
equilibrium phases can be treated with closely similar 
techniques. Nevertheless, this must be confirmed with test 
problems as in Paper I. 

2.8. A-iterations 

The robustness of the emissivities derived with the macro- 
atom formalism suggests that useful predictions of L v (t) 
may be obtained without converging to the exact NLTE 
solution. Instead, estimates of the excitation, ionization 
and temperature in each cell could be derived using the 
characteristics of the local MC radiation field as in Abbott 
& Lucy (1985) and in the previous diagnostic codes (Sect. 

I) . The accuracy of this approach in computing ioniza- 
tion fractions has been confirmed by Springmann & Puis 
(1998) for an O-star wind. 

Nevertheless, for definitive results, the NLTE solu- 
tion must be obtained. Fortunately, the use of indivisible 
^-packets and the macro-atom formalism facilitate this 
task. Because the thermal and statistical equilibrium con- 
straints are directly incorporated into the MC calculation, 
convergence to the NLTE solution can be achieved with 
geometry- independent A-iterations (Lucy 1999a; Paper 

II) . Thus, by simply repeating the process of bringing mat- 
ter into thermal and statistical equilibrium with the MC 
radiation field and then recomputing the latter, the solu- 
tion converges to the required equilibria. Moreover, con- 
vergence is rapid. 

With regard to thermal equilibrium, this success with 
A-iterations was initially demonstrated for a 1-D prob- 
lem (Lucy 1999a). Recently, the geometric independence 
of the technique has been demonstrated with its success- 
ful application to demanding 3-D problems (Harries et al. 
2004; Kurosawa et al. 2004). With regard to thermal and 
statistical equilibrium, experience is thus far limited to a 
simple 1-D problem (Paper II). 

2.9. Monte Carlo estimators 

A NLTE calculation requires the radiative rates of ex- 
citation, ionization and heating. With the RTE, these 
are obtained by numerical integration. But computing 
these quantities in a MC simulation is not so straight- 
forward. The obvious approach is simply to count the rel- 
evant events - e.g, photoionizations - in each cell in time 
At, thereby deriving the rate empirically. This is akin to 
a physics experiment in which detectors are distributed 
throughout an apparatus to record events. But, though 
appealing, this fails to make full use of the MC radiation 
field. For example, even if no photoionizations occur in a 
cell, a non-zero ionization rate must surely be derivable if 
it was traversed by r-packets containing ionizing photons. 

To derive MC estimators of radiative rates, a sum- 
mation procedure based on volume elements is preferred 



(Lucy 1999) to the more obvious choice of reference sur- 
faces. The basic building block from which the required 
estimators are derived will now be stated in a more gen- 
eral form than in earlier papers (Lucy 1999; Paper II). 

If, at a given position and time, I v (9, <j>) is the specific 
intensity at frequency v of a pencil of radiation propa- 
gating in direction (9,<p), then AttI v /c x dv duj is the in- 
stantaneous energy density of radiation in the frequency 
interval (y, v+dv) propagating within solid angle dco about 
the specified direction. This basic formula allows the MC 
radiation field to be converted into the conventional de- 
scription in terms of I v . 

If, during At, an r-packet of energy e„ with v € 
(v, v + dv) propagates in a cell of volume V, and if for a 
time St — Ss/c its trajectory is in the solid angle element 
du, then e^Ss/c is the packet's contribution to the cell's 
time- integrated radiant energy in dv dto. Accordingly, the 
volume- and time-averaged estimator of the specific inten- 
sity is given by 

^^ = r^E- fc ( 2 ) 

An AtV ^ e 

dfdui 

Here eo is the reference value for the energy of the £- 
packets (Sect. 3.2). 

Equation (2), when multiplied by the appropriate 
weight functions and integrated over v and lu, yields MC 
estimators for any required property of the radiation field. 
Estimators of this type were used initially to compute inte- 
grated mean intensities and absorption rates in a 1-D non- 
grey atmosphere in LTE (Lucy 1999a). Similar applica- 
tions to dusty circumstcllar envelopes have been reported 
by Wolf (2003) and by Niccolini et al. (2003). Estimators 
for the radiative rates needed for NLTE calculations can 
also be derived (Paper II). 

The class of MC estimators derived from Eq. (2) can be 
described as optimum and non-parametric. Optimum be- 
cause they use all the MC information and non-parametric 
because no assumption is made about the radiation field. 
These estimators are appropriate when the simulation is 
large enough that in each At every cell is traversed by 
many ^-packets. If not, a functional form can be assumed 
for I v or J v - e.g., a dilute black body - and its param- 
eters estimated from the limited number of trajectories 
in V. This functional form can then be multiplied by the 
appropriate weight functions to derive an estimate of the 
required quantity by numerical integration. This second 
procedure dampens the effects of sampling errors but, in- 
sofar as the functional form is inexact, yields biased es- 
timators - i.e., ones that do not converge to their exact 
values as Af — > oo, where M is the number of packets in 
the simulation. 

An application where the extra generality provided by 
Eq.(2) is essential is in constructing an estimator for the 
source function in a medium with non-isotropic scattering. 
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2.10. Emergent radiation 

The emergent rf spectrum for a MC simulation can be 
derived simply by counting escaping r-packets into fre- 
quency bins, with light travel-time taken into account as 
in Sect. 4.2. But even in the spherical case, large Af is then 
required to keep sampling errors small enough to allow a 
useful comparison with obseved spectra. The solution is to 
extract the source function from the simulation and then 
calculate the emergent flux from the formal integral (Lucy 
1999b). In a particular 1-D case, the errors in the resulting 
spectrum correspond to those of a binned spectrum from 
a simulation with J\f increased by a factor of <~ 320. For a 
3-D SN, where spectra for multiple lines-of-sight must be 
computed, the binning option is almost useless, and the 
already large gain factor with the formal integral will be 
vastly increased. 

3. Gamma ray transport 

In a previous paper (Paper II) on the NLTE transfer of 
UVOIR radiation in a SN envelope, the ^-packets were 
all created at the lower boundary, thus implicitly repre- 
senting the outward diffusion of the energy released by 
radioactive decays in the deeper interior. But with the 
computational domain now being the entire ejecta, the 
creation and transport of 7-rays must be treated explic- 
itly. In this section, therefore, the concept of indivisible 
^-packets is extended to cover this aspect of the general 
spectral synthesis problem for SNe. 

In a very recent paper, Milne et al. (2004) have com- 
pared results with various 7-ray transport codes and re- 
viewed the physical processes that must be treated. Their 
emphasis is on the prediction of 7-ray spectra rather on 
the powering of optical emission. 

3.1. Gamma ray lines 

The radioactive decays 56 Ni -> 56 Co and 56 Co -> 56 Fe 
each occur with the emission of a spectrum of 7-ray 
lines - see Table 1 in Ambwani & Sutherland (1988). 
Gamma-rays of energy Ei are emitted with probability 
/ ; when a parent nucleus decays. The total energy emit- 
ted per decay is therefore Elfu giving Eni = 1-728 and 
E Co = 3.566MeV for the 56 Ni and 56 Co nuclei, respec- 
tively. The e-folding times for these decays are ijvi = 8.80 
and tco = 113.7 days. 

The total 7-ray energy emitted by the decay sequence 
56 Ni -> 56 Co -> 56 Fe in the limit i -> 00 is 

Etot — {E Nl + Eco)M ra d/m Ni (3) 

where M ra d is the initial mass of 56 Ni, and mjvi is the 
mass of the 56 Ni nucleus. 

3.2. Radioactive pellets 

The initial mass of radioactive matter M ra d is quantized 
into N pellets, where N = E tot /eo. These pellets have the 



following property: they emit a single 7-packet with cmf 
energy eo and containing 7-rays whose cmf photon energy 
E corresponds to one of the lines emitted in the decay 
sequence 56 Ni — > 56 Co — > 56 Fe. Following this event, a 
pellet is permanently inert. 

As Af — > 00, this model must yield the correct time- 
dependent 7-ray line emissivities. To achieve this, we first 
identify two kinds of pellet: a fraction E^iKE^ + Ec ) 
are Ni pellets that collectively emit the 56 Ni spectrum, 
while the remainder are Co pellets, and they account for 
the 56 Co spectrum. 

If a pellet is designated as a Ni pellet, its decay time 
is randomly chosen as i 7 = —t^ln z; and the emitted 7- 
packet is assigned to line I from the 56 Ni spectrum with 
probability Eifi/J2 Eifi. On the other hand, a Co pellet's 
decay time is i 7 = —tNi^nzi — tc £nz2, where z\ and Z2 
are independent random numbers from (0,1). Two terms 
are required for the Co pellets since each 56 Co nucleus is 
created in the decay of 56 Ni nucleus. Having thus selected 
t 1 for a Co pellet, we select a line from the 56 Co spectrum 
as with 56 Ni. 

When a pellet decays, the emitted 7-packet is assigned 
a cmf direction vector fi in accordance with isotropic 
emission. Thus fi — (sin6 cos<fi, sin9 siruf>, cos9) , with 
cosO = 1 — 2z and <p — The corresponding rf vec- 
tor /i is obtained with the exact aberration formula - e.g., 
Castor (1972, Eq.(5)). The 7-packet's initial rf energy is 
then = eo/(l — n-v/c), where the local rf velocity 
v — r/t. 

The positions of the pellets at t\ are obtained by sam- 
pling the distribution of 5e Ni predicted by the explosion 
model. Their positions when they decay are then given by 
the assumption of homologous expansion. 

The 7-packets emitted in the interval (t n ,t n +i) are 
added to the 7- and r-packets still propagating in the 
ejecta at time t n to form a list of active ^-packets whose 
trajectories are up-dated during this time step. 

3.3. Transport of ^-packets 

In propagating through the ejecta, 7-rays create cascades 
of non-thermal electrons in multiple Compton scatterings 
off free and bound electrons as well as being absorbed in 
photoionizations. This redistribution of the energy of the 
emitted 7-ray into numerous channels over a substantial 
volume can nevertheless be modelled with indivisible 7- 
packets in such a way that the correct physics emerges as 
Af — > 00. 

3.3.1. Events 

As a 7-packet propagates, it undergoes events, both nu- 
merical and physical. The numerical events are: escaping 
from the grid, reaching the surface of a cell, or coming 
to the end of the current time step at t = t n+ \. The 
currently-included physical events are Compton scatter- 
ing and photoelectric absorption. 
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In describing how a 7-packet's trajectory is computed, 
it suffices to explain how to find the next event along the 
trajectory of one packet (cf. Paper II, Sect. 5). 

Given the rf position r and direction vector /x following 
an event, the next event is identified by computing the 
distances along the trajectory to all possible events and 
then selecting the event reached first. Since calculating 
distances to the numerical events is trivial, we here treat 
only the physical events. 

If E is the rf energy of photons in a 7-packet, the cmf 
energy is E = E (1 — fi.v/c). Accordingly, in the cmf the 
7-packet sees Compton scattering coefficient a (E ) and 
absorption coefficient k (E ). But in the rf these transform 
to a E = a (E )(l — fi.v/c) and k E — k (E )(1 — fi.v/c). 
Thus, when the radiative transport is carried out in the 
rf, the absorption, scattering (and emission) coefficients 
are direction-dependent - e.g., Castor (1972, Eqs. (2) and 
(3))- 

With these rf coefficients determined, we select the dis- 
tance 5s along the trajectory at which a physical event will 
occur from the standard MC formula 

(k E + <te) pSs = -tnz (4) 

and this event happens if 5s is smaller than the distances 
to the numerical events. Morever, when a physical event 
happens, it is a Compton scattering if z < o E j(k E + &e) 
and a photoelectic absorption if not. 

If 5s is the distance to the selected event, the updated 
space-time rf coordinates are r + 5s /x and t + 5s/c. 

3.3.2. Actions 

Starting with numerical events, we now describe subse- 
quent actions. 

If the 7-packet exits the grid, it escapes to 00, and 
attention then turns to the next 7-packet in the list of 
those active during the current time step. But if the 7- 
packet exits only the current cell and not the grid, then 
it enters the neighbouring cell and the search for the next 
event proceeds as above. Finally, if the event is the end 
of the time step at t n +i, computation of the trajectory 
is suspended, and the 7-packet's current rf data string 
r,t,[i,e E ,E is stored to await the next time step. 

With regard to physical events, the following actions 
are taken: In the event of a photoelectric absorption, the 
^-packet's trajectory as a 7-packet terminates. It now be- 
comes a k- or i-packet (cf. Sect. 2.7) with the same cmf 
energy e E as the absorbed 7-packet. Then, since energy 
storage in the gas is neglected (Sect. 2.4), this k- or i- 
packet converts immediately to an r-packet. In a general 
code, the frequency v of photons in this r-packet are de- 
termined by continuum emission (f-b or f-f) or line emis- 
sion following collisional excitation, and these are treat- 
able with the macro-atom formalism (Papers I & II). But 
here, with grey transport for UVOIR radiation, the emit- 
ted r-packet can be regarded as bolometric. Its cmf energy 
is e = e E and its cmf direction vector /x is selected ac- 
cording to isotropic emission. The rf direction vector /it 



then follows from the aberration formula, and the rf en- 
ergy is e = e /(l — n-v/c). These quantities together with 
r and t comprise the data string required to initiate the 
r-packet 's trajectory in the 3-D code for UVOIR radiation 
(Sect. 4). 

In the event of a Compton scattering, we first find 
the scattering angle 6 by randomly selecting cosQ from 
a look-up table of the percentiles of the cosQ probability 
distribution as functions of the incident 7-ray energy. A 7- 
ray scattered through angle 9 has its incident energy E = 
E m e c 2 reduced to fcE , where f c = l/[l+E (1— cos©)], 
with the remaining energy transferred to a Compton elec- 
ton. Thus the energy is divided in the ratio fc ■ 1 — fc- 
Accordingly, since ^-packets are indivisible, the 7-packct 
continues as a 7-packct if z < fc- If not, it becomes an 
e~ packet. 

If the packet continues as a 7-packet, its cmf energy re- 
mains e E but it now contains 7-rays of cmf energy fcE . 
To derive the rf values of these quantities, we must se- 
lect a new direction vector. In the cmf, the incident /x^ 
and emergent /x direction vectors are such that /x 1 ./x 
= cosQ. The combination of this with azimuthal angle 
$ = 2irz, where azimuth is defined with respect to polar 
direction , determines /x . The aberration formula then 
gives /x and so, in the rf, the emergent 7-packet has energy 
e E /(l — ijl.v/c) and contains 7-rays with photon energies 
fcE /(l — n-v/c). The next event for this 7-packet is now 
searched for as described in Sect. 3.3.1. 

On the other hand, if the 7-packct converts to an e~ 
packet, we assume in situ degradation into a bolometric r- 
packet with cmf energy e = e E . The rf data string needed 
to initiate its subsequent trajectory is now computed as 
for an r-packet created by photoelectric absorption - see 
above. 

3.4. Test 

The above treatment is a straightforward extension of the 
indivisible £-packet idea to 7-ray transport. Nevertheless, 
it should still be tested against the traditional MC treat- 
ment in which the MC quanta are photons (e.g. Colgate 
et al. 1980; Ambwani & Sutherland 1988). Such a test has 
been carried out by computing 7-ray deposition in a spher- 
ical SN. For simplicity, effects of 0{v/c) are neglected and 
photoelectric absorption acts only as a guillotine when E 
drops below lOOkeV. 

In one calculation, the pellets emit a single 7-packct 
as described above. In the comparison calculation, each 
pellet emits multiple photon packets with photon energies 
Ei and weights oc Ei fi , thus representing the appropriate 
7-ray line spectrum. Each such packet then undergoes sev- 
eral Compton scatterings, at each of which the energy of 
the Compton electron is deposited in situ. This continues 
until a photoelectric absorption deposits the remaining en- 
ergy. The two deposition profiles converge to one another 
as M — > 00. 
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3.5. Energy deposition rate 

In the present code, an explicit determination of H y , the 
heating rate due to the deposition of 7-ray energy is not 
required. The assumption of local thermal equilibrium im- 
plies that absorbed 7-packets are immediately re-emitted 
as r-packets (Sect. 3.3.2) and their subsequent propaga- 
tion is through matter with a temperature-independent 
grey absorption coefficient. But in a more general code, 
the thermal history of the ejecta must be calculated, and 
this requires 7i 7 . In any case, this quantity is needed here 
for the solution with moment equations (Sect. 5). 

In Sect. 3.3, two deposition processes are considered. 
First there is the loss of energy to Compton electrons, 
which occurs at the rate 

4 = 4tt J ' !{E)a{E)j E ,dE (5) 

where f(E) is the expected fraction of E transferred to 
the Compton electron and J E , is the cmf mean intensity 
at photon energy E . Second, there is the loss in photoion- 
izations, which occurs at the rate 

4 = 4tt j k{E)J E , dE (6) 

Estimators for these rates can be derived by applying 
Eq. (2) in the cmf. But note that the transport of packets 
is carried out in the rf. If 5s is the distance between con- 
secutive events for a 7-packet propagating in the rf, then, 
to 0(v/c), in the cmf, 5s = 5s (1 — fi.v/c). Accordingly, 
since we also have e E = ce (1 — H-v/c), the reqired esti- 
mators are 

c 'c = ^7 E £ -i 5s ^- 2/ ^/ c ) ( 7 ) 

and 

C 'i = S?7 E k ' ^ Ss (! - 2 ^/ c ) ( g ) 

These summations are over all trajectory segments 5s in 
V. Thus there is no restriction to packets that undergo 
Compton scatterings (Eq.[7]) or photoelectric absorptions 
(eq.[8]) in V. Accordingly, the resulting deposition rate 
7i 7 = C c + Cj is more accurate than the crude estimate 
obtained by summing the cmf energies of 7-packets that 
terminate their trajectories in V. 

3.6. Gamma-ray spectra 

A by-product of this 3-D time-dependent treatment of 7- 
ray transport is the prediction of 7-ray spectra, observa- 
tions of which have potential for detecting asymmetries 
in SN explosions (Hoflich 2002; Hungerford et al. 2003). 
Crude spectra can be obtained simply by binning escaping 
7-packets according to their rf energies E. But for high 
quality results, especially for the orientation-dependent 
spectra of an asymmetric SN, the formal integral should 
be used (Sect. 2.10). 



3.7. Future improvements 

Additional physical processes can be readily incorporated 
into this £-packet treatment of 7-ray transport. 

3.7.1. Positron transport 

Following Ambwani & Sutherland (1988, Table 1), 
positrons emitted in 56 Co are here assumed to annihilate 
in situ. But at late times, the transport of positrons should 
be followed until they escape or deposit their rest and ki- 
netic energies in the ejecta. This requires that we allow 
the Co pellets to emit e+ packets instead of the 7-packets 
with E = 0.511MeV. 

3.7.2. Transport of non-thermal electrons 

In treating Compton scattering, we assume that Compton 
electrons degrade in situ. The basis for this is the short 
stopping distance of MeV electrons compared to MeV 7- 
rays (Colgate et al. 1980). Nevertheless, we may eventually 
wish to follow the motion of the e~-packets, especially at 
late epochs. 

3.7.3. Pair production 

In addition to Compton scattering and photoelectric ab- 
sorption, a 7-ray with E > 2m e c 2 can tranform to an 
e + — e~ pair as it passes an atomic nucleus - see, e.g., 
Ambwani & Sutherland (1988) - with each of the pair 
having kinetic energy (ke) E/2 — m e c 2 . In stopping, the 
positron releases energy E/2+m e c 2 , comprising its ke plus 
the rest energy 2m e c 2 radiated when the positron annihi- 
lates with an ambient electron. In contrast, stopping the 
electron releases only its ke. 

Accordingly, when treated with indivisible ^-packets, 
pair production converts a 7-packet into an e + -packet with 
probability 1/2 + q or into an -packet with probability 
1/2 — q, where q = m e c 2 /E . In either case, the cmf energy 
of the packet is that of the incident 7-packet e E . If in situ 
deposition is not assumed, the selected packet's motion is 
then followed to escape or deposition. Note that a return 
to the pair production event to follow the other member 
of the pair is not required (cf. Paper I, Sect. 1). 

3.7.4. Rare radioactivites 

If the evolution of the ejecta is followed to late epochs 
(> 1000 days), then rare but long-lived radioactive iso- 
topes such as 57 Co, 44 Ti and 22 N a must be included (e.g., 
Woosley et al. 1989). This just requires additional types 
of pellet. 

4. Transport of UVOIR radiation 

In a previous paper (Paper II) the NLTE transfer of 
UVOIR radiation in a SN envelope of pure H was treated 
using and extending the macro-atom formalism (Paper I) . 
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Here, given the emphasis on time dependence and 3-D, we 
simplify to a grey absorption coefficient. 

4.1. Transport of r -packets 

As with 7-packets, the transport of r-packcts is carried 
out in a sequence of MC simulations for the time steps 
t n t n +i- Each of these simulations directly follows the 
corresponding one for 7-packets (Sect. 3.3). 

4.1.1. Sources of r-packets 

The r-packets that propagate in the ejecta in the time 
interval (t n , t n +\) have several sources. First are those that 
did not escape during the previous time step. For these, we 
have the rf data string (r,i,/z,e), and so their trajectories 
can be continued from time t n . Second are those created 
during the current time step as 7-packets are eliminated. 
This occurs either by photoelectric absorption or by the 
stopping of Compton electrons (Sect. 3.3.2). In either case, 
the 7 transport routine provides the rf data string (r,t,fi,e) 
needed to initiate their trajectories as r-packets. 

A third source is specific to the first time step. Radiant 
energy emitted by pellets that decay at times prior to ti 
must be accounted for. An approximate treatment is as 
follows: A 7-packct emitted at i 7 < ti and position r 7 
is assumed to have converted to an r-packet by time t\ 
but to have diffused negligibly relative to matter (posi- 
tion coupling) in the time interval (t 7 ,£i) so that at t\ 
its rf position is r 7 (fi/i 7 ). This is justified by the short 
mean free paths of photons in these early, high-density 
phases. The rf direction vector at t\ is computed on the 
assumption of isotropic emission in the cmf. 

The energies of these packets must also be specified. 
Although position coupled, they still do work on the ex- 
panding ejecta. An r-packet's energy at t\ is therefore 
e = €q tj/ti. 

An additional source of r-packets at t\ is the radia- 
tion generated by shock heating during the explosion. The 
explosion model's prediction of this radiation can be dis- 
cretized into r-packets of energy e 0j the same quantum of 
energy as for the radioactive pellets (Sect. 3.2). However, 
such r-packets are neglected in this test code. 

4.1.2. Events 

As an r-packet propagates, it undergoes numerical and 
physical events. The numerical events are identical to 
those for 7-packets (Sect. 3.3.1). The physical events are 
absorptions. 

In describing how a r-packet's trajectory is computed, 
it suffices to explain how to find the next event along the 
trajectory of one packet. Given the rf data string (r,t,n,e), 
the next event is identified by computing the distances 
along the trajectory to all possible events and then select- 
ing the event reached first. As with 7-packets, we treat 
only physical events. 



Because the transport of r-packets is also carried out in 
the rf, the anisotropy of the rf absorption coefficient must 
again be allowed for, as it was for 7-packets in Sect. 3. 3.1. 
If the grey absorption coefficient per unit volume in the 
cmf is k , the effective value in the rf is k = k (1 — /z.u/c), 
and so the distance 5s to the absorption event is given 
by k Ss = —inz. This event happens if this is smaller 
than the distance to any numerical event. With Ss thus 
determined, the coordinates (r, t) are updated as for 7- 
packets (Sect. 3.3.1). 

4.1.3. Actions 

For numerical events, the subsequent actions correspond 
to those for 7-packets (Sect. 3.3.2). 

If the event is an absorption, we assume instantaneous 
isotropic re-emission in the cmf. Accordingly, the new rf 
direction vector /x is calculated as it was for emitted 7- 
packets in Sect. 3.2. 

The energy content of the bolometric r-packet must 
also be updated. If the incident packet has rf energy e\ and 
direction vector fii, its cmf energy is e = ei(l — fi^.v/c). 
Since this is conserved by the absorption-emission event, 
the updated rf energy is e = e /(l — fi.v/c). Calculation of 
the post-event data string (r,i,/x,e) is now complete, and 
the search for the next event can begin. 

4.2. Bolometric light curve 

Consider an r-packet with post-event data string (r,t,/x,e) 
that is an escapee from the grid. A distant observer at rest 
in the rf who detects this packet records its arrival at what 
he perceives to be a time r = t — fi.r/c after the explosion. 
Thus the pairs (e, t) detected by one such observer is the 
data set from which he can construct the bolometric light 
curve of the 3-D SN as seen from his orientation. 

In this paper, the 3-D code is tested by applying it to 
a spherically-symmetric SN. Accordingly, we use all pairs 
(e, t) to construct its orientation-averaged bolometric light 
curve (Sect. 6.3). 

5. Solution with moment equations 

In order to test the 3-D code described in Sects. 3 and 
4, we apply it in Sect. 6 to a spherically-symmetric SN. 
This allows the code to be tested against a solution of the 
same problem obtained by numerical integration of the 
time-dependent RTE. 

5.1. Castor's equations 

Castor (1972) has given a general treatment, accurate 
to 0(v/c), of radiative transfer in spherically-symmetric 
flows. In particular, he derives the zeroth and first 
frequency-integrated moment equations in the cmf. These 
two equations, applied to a homologously expanding flow 
with grey absorption, are the basis of the comparison cal- 
culation. 
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If we again neglect the contributions of the gas to en- 
ergy balance (Sect. 2.4), the zeroth moment equation is 



dt 



+ P 



dL 



(9) 



where the dependent variables are U R , the cmf energy 
density of radiation, and L , the cmf luminosity variable. 
The independent variables are elapsed time t and the mass 
coordinate M r . Note that time derivatives in the moment 
equations are Lagrangian. 

The first moment equation is 



1 dL dP R 

c dt oM, 



+ a 2 (3P R -U R ) = -(k +-)L (10) 



dt' 



where a\ = 167r 2 r 4 pc, a 2 = 4nrc, and the third dependent 
variable P R is the cmf radiation pressure. 

5.2. Closure approximation 

To make this system determinate, an approximate formula 
relating the three moments must be imposed. Here we 
adopt Eddington's approximation, 



Pn 



which is used to eliminate P R from Eq.(10). This leaves 
two partial differential equations (PDEs) in the two vari- 
ables U R (M r ,t) and L R (M r ,t). 

5.3. Boundary conditions 

The PDEs must be solved subject to appropriate bound- 
ary conditions. Clearly, at the SN's centre, 



L (0,t) = 



(12) 



while, at the surface, Eddington's boundary condition 
F(0) = 2J(0) gives 



L (M,t) = 2nR 2 max cU R (M,t) 



(13) 



Here R m ax(t) = v max t, and M is the total mass of the 
eject a. 

5.4. Initial conditions 

In addition, initial conditions are required at t = t\. As 
in the MC code (Sect. 4.1.1), we assume that diffusion of 
radiation relative to matter is still negligible (Sect. 4.1.1) 
at ti, and so 



L (M r ,h) = 



(14) 



An initial condition for U R is derived as follows: For 
t < t\ , negligible relative diffusion of 7-rays implies in situ 
deposition, so that 7i 7 = 47rj 7 , the integrated 7-ray emis- 
sivity in the cmf. Accordingly, from Eq.(9), the required 
initial condition is given by integrating the ordinary dif- 
ferential equation (ODE) 



dUz 
dt 



+ 



4U 



R _ 



= 4tt1 



which can be done analytically. 

If f{M r ) is the mass fraction of 56 Ni in the mass shell 
M r at t = 0, then 



4nj = S(t)f p/m m 



(16) 



where S(t), the energy released per 56 Ni nucleus, is given 

by 



_ Em t/t Ni 



S(t) = —e 
tm 



E C o 



tCo — tw. 



(e 



-t/t c 



-t/tNi 



) (17) 



Since each term in Eq. (17) has the same form and 
Eq. (15) is linear, we consider a single exponential - i.e., 
S(t) = E t /Ue-*' u , and we note also that homologous 
expansion implies that p cx 1/t 3 . Integration of Eq.(15) 
then gives 



E* 



t 



[1 



-t/t, 



t 



f pj m Nl 



(18) 



This solution can be applied to each term in Eq. (17) to 
construct the complete solution. But with t\ <§; tm, the 
first ( 56 Ni) term suffices. 

Note that, in deriving Eq.(18), we assumed 



(11) U R (M r ,t) = at t 



0. This is consistent with the 
neglect of shock heating in the MC code (Sect. 4.1.1). 

5.5. Solution technique 

Eqs. (9) and (10) are solved in the same way as are the 
equations of stellar evolution when the P and T terms are 
included in the energy equation (e.g., Schwarzschild 1958, 
p. 100). Thus, the time derivatives of U R and L are re- 
placed by backward difference formulae. Then, since the 
solution at earlier time steps is known, the PDE problem is 
effectively simplified to a two-point boundary value prob- 
lem for ODEs. When the space derivatives in the ODEs 
are approximated by difference formulae, the resulting al- 
gebraic system can be solved with an elimination proce- 
dure (Henyey method). Here, a slight generalization of the 
procedure described by Richtmyer (1957, pp. 101-104) has 
been followed. 

In stellar evolution codes, it is still current practice (A. 
Weiss, private communication) to approximate the time 
derivative of a variable Q at time t n with the formula 



dQ 
dt 



AQ r , 



At n . 



(19) 



where A denotes the forward difference operator, so that 
Ai n _i = t n — t n -\. This formula is not centred and thus 
has error O(At). But higher accuracy can be achieved 
by using more than one previous model - c.g, Richtmyer 
(1957, p. 94, item 9). If the two previous models at t n -\ 
and t n -2 are used, the appropriate difference approxima- 
tion is obtained by fitting a parobola and then evaluating 
its derivative at t n . The resulting formula is 



At„_ 



At n - 



(20) 
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where a — At„_i/(At„_i + At n -<i). Thus, with constant 
time steps, a = 1/2. The gain in accuracy is investigated 
in Sect. 6.2. 

In starting the integration of the PDEs, the theory of 
Sect. 5.4 is used to provide the previous solutions required 
by Eqs. (19) and (20) as well as initial estimates of U R and 
L for the first Henyey iteration. 

5.6. Light curve 

The solution of the PDEs yields the two functions 
U R {M r ,t) and L (M r ,t). Thus we directly get the bolo- 
metric light curve L (t) seen by a cmf observer at A4 r = 
A4, the surface of the SN. But the quantity of interest is 
the light curve seen by a distant observer at rest in the rf. 
This can be calculated using ^-packets as follows: 

In a small interval At at t, the surface emits N r- 
packets at random times t n in At. These all have cmf en- 
ergy e = L (t)At/N, and their direction cosines are jj, n = 
\fz, corresponding to zero limb-darkening, an assumption 
consistent with Eddington's approximations. Their direc- 
tion cosines fi n in the rf are given by the aberration for- 
mula, and their rf energies are then e„ = e/(l — ^ n v ma x / c) . 
The distant observer perceives the nth packet as having 
been emitted at elapsed time t„ = t n — PnRmax/c, and so 
the N pairs (e„, r„) represent the contribution of the inter- 
val At to the bolometric light curve seen by this observer. 
Summing over all intervals At then gives the complete 
light curve (cf. Sect. 4.2). 

6. Numerical results 

After first investigating the accuracy of bolometric light 
curves for spherical SN obtained with the cmf moment 
equations, this section uses such a light curve to check the 
3-D MC code described in Sects. 3 and 4. 

6.1. Model 

The model SN is closely similar to that used by Pinto 
& Eastman (2000) to investigate the parameter sensi- 
tivity of type la light curves. The ejecta has uniform 
density and its basic parameters are: M. = 1.39.Mq, 
M( 56 Ni) = 0.625Mq, and v max = 10 4 km s" 1 . The 
56 Ni is assumed to be stongly concentrated in the cen- 
tral core. Thus, f(M r ) = 1 for M r < O.5A^0 and then 
drops linearly to zero at A4 r — 0.75Mq. 

The grey absorption coefficient for UVOIR radiation is 
k I p — 0.1 cm 2 g _1 , and the photoelectric absorption co- 
efficient k E for 7-rays is derived from Eq. (3) of Ambwani 
& Sutherland (1988) with nuclear charge Z = 14. 

6.2. Accuracy of moment solution 

Numerical solutions of Castor's equations have been ob- 
tained in order to test the MC code against an RTE treat- 
ment accurate to O(v/o). But the neglect of higher order 
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Fig. 1. Bolometric magnitude at maximum as function of 
time step A logt. Time derivatives are approximated with 
Eq. (19) {open circles) or Eq.(20) {filled circles). The ex- 
trapolation to vanishingly small time step is shown {star) 
as are ± 0.05 mag. error bands. 

terms in v/c is not the only source of error in the RTE 
solutions. Larger errors may arise due to the closure ap- 
proximation (Sect. 5.2) and the difference approximation 
of time derivatives (Sect. 5.5). 

The moment equations are solved on a uniformly- 
spaced grid with 400 grid points, the innermost at 
10~ 3 R max . The energy deposition rate H 7 in Eq. (9) is 
derived with a 1-D version of the MC code described in 
Sect. 3 using the estimators given in Eqs. (7) and (8). The 
number of radioactive pellets is J\f = 10 7 . 

To investigate sensitivity to the approximation of time 
derivatives, two sequences of light curves are computed, 
one using Eq.(19) to evaluate dU R /dt and dL R /dt, the 
other using Eq.(20). Along each sequence, A logt varies, 
thereby determining the rate of convergence as A log t — > 
0. For each light curve, M^ x , the peak rf bolometric mag- 
nitude, is obtained by parabolic fitting. 

The values of M^ x as a function of A logt are plotted 
in Fig. 1, together with the value M™f x = -19.181 ob- 
tained by extrapolating to infinitesimal time step. Error 
bands of ±0.05 mag. are also drawn. 

Fig. 1 shows, as expected, that errors grow linearly 
with Eq. (19) and quadratically with Eq. (20). From 
this plot, we find that accuracy < 0.05mag. requires 
A logt < 0.048 with Eq.(20) but A logt < 0.009 with 
Eq.(19). Clearly, for high precision, Eq (20) is preferred 
and is indeed used in Sect. 6.3. 

Although this test refers to the cmf moment equations, 
it is surely relevant generally to numerical solutions of the 
time-dependent RTE. Accurate difference approximations 
for time derivatives are necessary to limit the accumula- 
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Fig. 2. Comparison of bolometric light curves. The light 
curve obtained with the 3-D Monte Carlo code (filled cir- 
cles) compared to that obtained from the moment equa- 
tions (solid line). The corresponding 7-ray deposition 
curves are also plotted. 
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Fig. 3. Light curve residuals. The differences between the 
bolometric light curves plotted in Fig. 2 as a function of 
elapsed time t. Residuals for 3-D MC calculations with 
100 3 (filled circles) and 20 3 (open circles) grids are shown. 



tion of errors as one integrates forward in time. And of 
course the same remark applies to stellar evolution codes. 

6.3. Bolometric light curves 

In Fig. 2, the rf bolometric light curve computed with the 
3-D MC code is compared with that derived from the cmf 
moment equations. The MC light curve is from a simula- 
tion in which the initial distribution of radioactive matter 
is represented by TV = 4 x 10 6 pellets, with the resulting 
7- and r-packets propagating in a 100 3 grid. The calcu- 
lation starts at logt\(dys) = 0.3, and the time steps are 
A logt — 0.01. The rf light curve is derived from escaping 
r-packets as described in Sect. 4.2. Thus, for packets with 
arrival times r in the interval (t n ,t n+ i), the values of e 
are summed to obtain an estimate of L n+ i/ 2 . 

The rf light curve derived from the PDEs of Sect. 5.1 
is also plotted in Fig. 2. This is derived from the cmf light 
curve L (t) as described in Sect. 5.6. As for the calculations 
of Sect. 6.2, the spherical SN is modelled with 400 shells, 
and Af = 10 7 in the 1-D 7-ray code. 

Fig. 2 shows that the two bolometric light curves are 
in good agreement over an extended time interval, from 
well before to long after maximum light at t = 15.3 days. 
This implies a corresponding degree of agreement in the 
7-ray energy deposition rates. Nevertheless, these are also 
plotted in Fig. 2. For the 3-D code, the deposition rate is 
calculated simply by summing the energies of 7-packets 
as they convert to r-packets. But for the 1-D code, the 
estimators defined by Eqs. (7) and (8) are used to calculate 
H , which is then summed over shells. 



To investigate the precision of the 3-D light curve in 
more detail, the residuals relative to the PDE solution 
are plotted in Fig. 3. This plot reveals large systematic 
residuals at early times on the rising branch of the light 
curve but that these decrease with time becoming quite 
small for t <; 8dys. For the interval 10 — 50 days, the mean 
residual is -0.008 mag. A residual of this order may reflect 
errors in the PDE solution since this relies on Eddington's 
approximations. 

The large residuals on the rising branch have been 
traced to inadequate spatial resolution of the distribution 
of 56 Ni in the 3-D code, despite the 100 3 grid. Because 
56 Ni is uniformly distributed within each cubical cell, the 
function f(Ai r ) (Sect. 6.1) is slightly broadened relative 
to its more accurate representation in the 1-D calculation. 
This results in some released radioactive energy reaching 
the surface slightly early, giving a brighter MC light curve. 
This explanation is confirmed by finding that these resid- 
uals become larger with a coarser grid. To illustrate this, 
Fig. 3 includes the residuals for t < 10 days with a 20 3 
grid. 

This comparison with the PDE solution shows that, 
with a fine enough grid, the MC code described in Sects. 
3 and 4 is capable of carrying out high precision transport 
calculations accurate to 0(v/c) for 7-rays and UVOIR 
radiation propagating in a 3-D SN. 

6.4. Diagnostics 

Among the merits of using ^-packets is conservation of en- 
ergy to high precision and the ready monitoring of energy 
transformations within the configuration. 
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Fig. 4. Energy fractions as functions of elapsed time t. 
Quantities plotted are E 1 , the integrated energy released 
by radioactive decays, E^, the integrated energy emitted 
to oo, W, the work done on the expanding ejecta, and En, 
the radiant energy stored in the ejecta. The unit of energy 
is E to t - see Eq. (3) - so that E 7 (oo) = 1 . 

6.4.1. Energy conservation 

Energy conservation for the ejecta implies that 

E 00 (t) + E R (t)-W(t)=E 1 (t) (21) 

Here E 7 (t) is the total energy released by radioactive de- 
cays in the time interval (0,t), £ , 00 (i) is the total energy 
lost through the surface in (0,t), and E R (t) is the radiant 
energy stored in the ejecta at time t, all these quantities 
being evaluated in the rf. Finally, W(t) is the total work 
done in (0, t) by radiation interacting with the expanding 
ejecta. 

Estimators for the quantities in Eq. (21) are sums over 
packets. Thus, E 1 (t) is the sum of the initial rf energies 
ee of 7-packets emitted by radioactive pellets (Sect. 3.2); 
Eoo(t) is the sum of the rf energies eg and e of 7- and r- 
packets that escape the grid in (0,t); and En(t) is the sum 
of the rf energies of packets still propagating within the 
ejecta at time t. Finally, W(t) is the sum over all physical 
events in (0,t) of Ae = e\ — 62, where ei and €2 are an 
^-packet's incident and emergent rf energies, respectively. 

At all time steps in the calculation of Sect. 6.3, the 
left- and right-hand sides of Eq. (21) agree to better than 
1 part in 10 12 . 

6.4.2. Energy flows 

The time variations of the quantities in Eq. (21) are plot- 
ted in Fig. 4, where the unit of energy is E tot from Eq. (3). 
For the adopted parameters, E to t = 1-14 x 10 50 erg. 



Fig. 4 illustrates the strong departures from station- 
arity. In the early dense phases (t ^ 5 days), energy re- 
leased by radioactive decays is trapped within the ejecta, 
resulting in increasing E R (t) while E^t) remains close 
to zero. This continues until En(t) reaches a maximum of 
0.085 at t = 9.6 days. Thereafter, the dropping density 
allows the trapped radiation to be released, resulting in 
a sharp increase in E^t) and concomitant decrease in 
Eji{t). For t <; 35 days, energy storage is inconsequential 
(En ^ 0.02); stationarity then becomes a good approx- 
imation, as is evident from the closely similar slopes of 
E 7 (t) and 

The maximum of En (t) is a rough measure of the high- 
est fraction of ^-packets active at one time and for which, 
therefore, storage space must be allocated. In fact, be- 
cause many of the packets trapped at t = 9.6 days have 
lost energy by doing work on the expanding ejecta, the 
active fraction is larger. In the calculation reported here, 
the peak active number is 7.2 x 10 5 at t = 9.2 days - i.e., 
0.18JV. 

Non-stationarity is also illustrated by Fig. 5, which 
plots residence against escape times for ^-packets in a 
small simulation with TV = 2000. For t & 20 days, we see 
that it is not uncommon for escaping r-packets to have 
residence times close to the age of the SN. Thereafter, the 
lowered density allows packets to escape readily and typ- 
ical residence times drop to ~ a few days. Note also that 
no 7-packets escape until t 20 days. 

In stationary problems, or this problem with c = 00, 
a MC simulation must continue until all packets have es- 
caped. This limits the usefulness of MC methods for op- 
tically thick media. But here a packet initially at large 
optical depth escapes when that depth has dropped to 
~ 1 because of expansion. Fig. 5 illustrates this effect. 

For this problem, the difficulty with an unacceptably 
large number of events before escape would reappear if 
t\ — > were necessary for accuracy. Fortunately, the posi- 
tion coupling assumption (Sect. 4.1.1) is well justified for 
t ^ 1 day. 

7. Conclusion 

The aim of this paper has been to initiate a MC ap- 
proach to computing light curves and spectra for 3-D 
SN explosions. The adopted procedure is based on in- 
destructible and indivisible ^-packets as the MC quanta. 
Although not called upon for this test problem, the at- 
traction of this approach in the long-term for this and 
other multi-dimensional NLTE problems is that it al- 
lows the constraints of thermal and statistical equilibrium 
to be incorporated into the scheme via the concept of 
macro-atoms. Moreover, the iterations required to achieve 
full self-consistency between radiation, level populations 
and the electron temperature can be carried out with 
geometry- independent A- iterations . 

The calculations reported in Sect. 6 demonstrate that 
high photometric precision can be achieved with MC 
methods for time-dependent transport calculations in 3- 
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each time step for full self-consistency should initially be 
omitted in favour of approximate formulae relating gas 
characteristics to the local MC radiation field. 



Acknowledgements. I am grateful to S.A.Sim for detailed com- 
ments on the manuscript. 



Fig. 5. Residence times At — t esc — t 1 plotted against 
t esc . Escaping r packets ( open circles) and 7-packets {filled 
circles) are indicated. 



D. Contributing strongly to maintaining accuracy are the 
radioactive pellets introduced in Sect. 3.2 since they result 
in a seamless transition from energy transport by 7-rays 
to that by UVOIR radiation. This avoids the loss of accu- 
racy that might arise if the 7-ray deposition profile were 
separately calculated and then randomly sampled to cre- 
ate UVOIR radiation. An incidental benefit is simplified 
coding. 

In Sect. 1, this investigation was described as provid- 
ing a platform onto which more detailed treatments of 
matter-radiation interactions can be added. Thus the al- 
ready detailed modelling of 7-ray transport can be further 
improved as discussed in Sect. 3.7. Also C c , the rate of 
energy deposition in the form of Compton electrons, is 
available from Eq. (7) as the source term for a Spencer- 
Fano calculation of the rates at which non-thermal elec- 
trons lose energy to thermal electrons and in ionizing and 
exciting atoms. Finally, the NLTE transfer of UVOIR ra- 
diation can be implemented following the successful test 
for a pure H envelope in Paper II. 

An encouraging aspect of this agenda is that the 
physics missing from the present code is well understood, 
even though many relevant cross sections are still poorly 
known. But a discouraging aspect is that the demands 
on computer power are such that useful results for a 
physically- realistic 3-D model are not feasible at present 
with a single workstation. Nevertheless, as argued in Paper 
II, an attractive option is to implement these methods on 
a computer with numerous parallel processors. But even 
then, it is probably desirable and necessary to approach 
the full problem via simplified versions that are less de- 
manding on computer time and storage space. Specifically, 
as in previous diagnostic codes, the iterations required at 
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